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Abstract 



Variational data assimilation technique applied to identification of optimal approximations 
OO ■ of derivatives near boundary is discussed in frames of one-dimensional wave equation. 

Simplicity of the equation and of its numerical scheme allows us to discuss in detail as 
the development of the adjoint model and assimilation results. It is shown what kind of 
errors can be corrected by this control and how these errors are corrected. This study is 
carried out in view of using this control to identify optimal numerical schemes in coastal 
regions of ocean models. 



Key words: Variational Data Assimilation, Boundary conditions; Wave equation 
PACS: 47.85.L 

. 1. Introduction. 

^ . It is now well known, even the best model is not sufficient to make a good forecast. 

\ Any model depends on a number of parameters, requires initial and boundary conditions 

O ' and other data that must be collected and used in the model. However, interpolating or 
0\ ' 
O 



smoothing observed data is not the best way to incorporate these data in a model. Lorenz, 
in his pioneer work [l| has shown that a geophysical fluid is extremely sensitive to initial 
conditions. This fact requires to bring the model and its initial data together, in order to 
work with the couple "model-data" and to identify the optimal initial data for the model 
5^ \ taking into account simultaneously the information contained in the observational data 
and in the equations of the model. 

Optimal control methods [2| and perturbations theory [3] applied to the data assimi- 
lation technique ([^, j^) show the way to do it. They allow to retrieve an optimal initial 
point for a given model from heterogeneous observation fields. Since the early 1990 's, 
many mathematical and geophysical teams have been involved in the development of the 
data assimilation strategy. One can cite many papers devoted to this problem, as in the 
domain of development of different techniques for the data assimilation and in the domain 
of its applications to the atmosphere and oceans. 

However, overwhelming majority of data assimilation methods are now intended to 
identify and reconstruct an optimal initial point for the model. Since Lorenz who 
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has pointed out the importance of precise knowledge of the starting point of the model, 
essentially the starting point is considered as the control parameter and the target of the 
data assimilation. 

Of course, the model's flow is extremely sensitive to its initial point. But, it is rea- 
sonable to suppose that a geophysical model is also sensitive to many other parameters 
like bottom topography, boundary conditions on rigid and open boundaries, forcing fields 
and friction coefficients. All these parameters and values are also extracted in some way 
from observational data, interpolated to the model's grid and can neither be considered as 
exact, nor as optimal to the model. On the other hand, due to non-linearity and intrinsic 
instability of model's trajectory, its sensitivity to all these external parameters may also 
be exponential. 

Numerous studies show strong dependence of the model's flow on the boundary data 
on the representation of the bottom topography (jsl, j^, 10]), on the wind 
stress ([lH, l3]), on diffusivity coefficients ( l3) and on fundamental parametrization 
like Boussinesq and hydrostatic hypotheses [IJ]. But few papers are devoted to the 
development of data assimilation techniques intended to identify and to control these 
model's parameters. One can cite several attempts to use data assimilation in order to 



identify the bottom topography of simple models f|15 |. jl6|] ') and in order to control open 



boundary conditions in coastal and regional models ([17[, [l8j, [19|). Boundary conditions 
on rigid boundaries have been controlled by data assimilaton for heat equation (see for 
example 2^], 21]), but this control concerns the diffusion operator rather than transport 
and advection type operators used in geophysical models. 

This paper presents a preliminary study of using variational data assimilation in order 
to identify an optimal parametrization of boundary flows and boundary conditions on 
rigid boundaries. Despite the boundary configuration of the ocean is steady and can 
be measured with much better accuracy than the model's initial state, it is not obvious 
how to represent it on the model's grid because of limited resolution. The coastal line of 
continents possesses a very fine structure and can only be roughly approximated by the 
model's grid. Consequently, boundary conditions are defined at the model grid's points 
which are different from the coast. Even the most evident impermeability condition being 
placed at a wrong point may lead to some error in the model's solution. From physical 
point of view, we should accept the fiux can cross the boundary in places where the 
boundary is in water, prescribing some integral properties on the fiux. 

Even in a fine resolution model when boundary currents are explicitly resolved, it is 
not clear what kind of boundary conditions to prescribe for tangential velocities. However, 
prescribing slip or no-slip conditions may result in a drastic change of the global circulation 
(see (sl). 

Consequently, it may be reasonable to use variational data assimilation in order to 
determine what boundary conditions are optimal for the model's variables. However, 
instead of controlling boundary conditions themself, it may be more useful to identify 
optimal discretization of differential operators in points adjacent to boundaries because 
this is more general case. Indeed, boundary conditions participate in discretized operators, 
but considering the discretization itself, we take into account additional parameters like 
the position of the boundary, lack of resolution of the grid, etc. 

In this paper we use data assimilation to control the discretization of derivatives in 
adjacent to boundary grid points. The development of the data assimilation is illustrated 
on the example of the simplest one-dimensional wave equation. On one hand, the sim- 
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plicity of the equation allows us to clearly see technical points of the development (like 
the algorithm of differentiation and development of the adjoint equation) without being 
overwhelmed by complexity of operators and grids. On the other hand, the knowledge 
of the exact solution and of the errors of numerical discretization of the wave equation 
allow us to clearly see how these errors are corrected by data assimilation. The purpose 
of the paper is to study the possibility to control boundary numerical scheme by data 
assimilation and the particularities of this type of control in view to develop and use the 
data assimilation to identify optimal numerical scheme in coastal regions of ocean models. 

The paper is organized as follows. The second section describes the model, its ad- 
joint and the data assimilation procedure. The third section is devoted to numerical 
experiments and discussion. 

2. One-dimensional wave equation 

As it has been noted in the introduction, we consider one-dimensional wave equation 
written for u — u{x, t) and p — p{x, t) in the following way: 

du dp Q 

dt dx 

I - I- w 

This equation is defined on the interval < x < 1 with boundary conditions prescribed 
for u only: 

u{0,t) ^u{l,t) ^0 (2) 
Initial conditions are prescribed for both u and p 

u{x, 0) = u, p{x, 0) = p (3) 

The equation is discretized on a regular grid that is somewhat similar to Arakawa's C 
grid in two dimensions: 

Ui = u{ih) for i = 0, . . . N 
Pi_i/2 = p{ih -h/2) for t = l,...N (4) 

with h = This grid is well adapted to the prescribed boundary conditions because 
the boundary points x — and x — 1 belong to the grid for u discretization, but do not 
belong to p-grid. 



Pl/2 P3/2 P5/2 PN-5/2 PN-3/2 PN~1/2 

I ^ 1 ^ 1 X 1 1 X 1 X 1 X — I 

Uo Ui U2 U3 Un-3 Un-2 UN-1 Un 



Discrete derivatives of u and p are defined as follows 



at all internal points in the interval, i.e. 2 < i < N ~ 2 for and 1 < i < — 2 

for • Coefficients aj are supposed to be known because we intend to control 

approximations near the boundary only. In this paper, we use either the sequence aj = 
(0, —1, 1, 0) or the sequence aj = ^(1, —27, 27, —1) for j = (—1, 0, 1, 2). One can easily 
see that corresponding approximations are of second and of fourth order approximation 



h 



p,_3/2 - 27pi_i/2 + 27pi+i/2 - Pi+3/2 ^ fdp\ _3h^/d^\ 



(6) 



To be able to solve numerically the equation ([T]), we need also to approximate deriva- 
tives of u and p near boundaries at points i = 1/2, N — 1/2 and i = 1, N — 1 respectively. 
These approximations are supposed to be different from (jS]) and include the control vari- 
ables in this study. Moreover, expressions can not be used at all for the fourth order 
approximation because they require function's values beyond the boundary: u^i and 
P-i/2- We can, of course, extrapolate u and p beyond the domain with the necessary 
order and substitute extrapolated values in ([5]), but it is not obvious what extrapolation 
formula is the best for this purpose, especially for p. So, in order to obtain an optimal 
boundary approximation assimilating external data, we suppose nothing about derivatives 
near the boundary points and write them in a general form 

dp\ l>r^ p 

We do not fix the value of J in these formula intentionally because we shall see further 
its infiuence. 

Here we can emphasize the choice to control the numerical scheme in the boundary 
region rather than boundary conditions. The general form of boundary conditions that 
may be prescribed for u variable of the one dimensional wave equation writes 

u{0,t)-A—{0,t)=B. 

We can not impose more complex boundary conditions (with second derivatives, for ex- 
ample) because we obtain a system with no solution at all. Consequently, we can control 
only two parameters, A and B. It may be sufficient in particular cases, but, as we shall see 
further, is not sufficient in general. However, controlling all coefficients of the numerical 
scheme ([7]) , we are free to choose as many aj as we need defining appropriate value of the 
parameter J. 

We distinguish and a" allowing different derivatives approximations for p and for 
u because of the different nature of these two functions and different boundary conditions 
prescribed for them. Derivatives at the opposite side are calculated by 

dp\ 1 ^ 

j=o 



dx/N-i h 
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dxj 



N~l/2 
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and coefficients and a" are also considered as unknown control parameters different 
from and a". All together, we have 4( J + 1) control parameters. 
Time stepping is performed by leap-frog scheme 



2t 



dxJi 



„n+l _ „n-l 
Pi-1/2 Pi-1/2 

27 



dx) i-l/2 







(9) 



The ffist time step is splitted into two stages in order to ensure second order approximation 
in time and to avoid typical leap-frog splitting between odd and even timesteps. 



1/2 



r/2 



dx) i 



u]-u^ 



dx 



^ =0, 



1/2 
Pi-1/2 ~ Pi-1/2 

?72 



ox/ i-l/2 



Pi-1/2 ^ P^-1/2 



dx J i-l/2 



0. 
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Approximation of the derivative introduced by ([5]) and ([7j) depends on control variables 
a. The operator is not completely defined as in usual schemes, but it is allowed to change 
its properties near boundaries in order to find the best fit with requirements of the model 
and data. To assign variables and a" we shall perform data assimilation procedure 
and find their optimal values. 



2.1. Tangent and adjoint equations 

First of all, we calculate the Gateaux derivative of the operator with respect to control 
parameters. Control variables are supposed to have small variations and we determine 
how these variations will perturb the solution of the model. Thus, we suppose that all a 
are replaced by some a + Sa such that \\Sa\\ << ||a||. Let the model with a + Sa have a 
new solution u + 6u, p + 6p. In this case, variables 6u, Sp must satisfy 



d6u 
dSp 



- D^P'^Sp - SD^P'^p - SD^P'^Sp = 



where operators /^^^^(a^) and D^^\a^) are approximations of derivatives defined by 
([7]), and ([8]), i.e. for the p derivative 
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Operators 5D^^^ and 5D'^'^^ are the differences 
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(13) 



and similarly for operators and D^^\ 

However, expressions 6D^^^p and 5D^'^\ in ( ITTl) are not convenient to make further 
development. Writing an adjoint operator, we would better have a constant operator, 
which does not depend on 6a, multiplied by a variable vector which depends on 6a. This 
is the case in products D^^6p where 6p depends on 6a, but it is not the case in products 
like 6D^P^p where p is solution of original equation and has no relation with 6a. It would 
be more convenient to rewrite these products: 



/ 



6D^P^ 



p 



1 
h 



j=0 



E 6a^jPj+i/2 






V 



j=0 



E 6a^jPN- 



-i-i/2 



P6aP 
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3=0 







' E 6a'^UN~ 

j=0 
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U6a'' 



(14) 



where operators P and U are constructed from the solution p and u of the original 
equation. Their matrices have non-zero elements in the first and in the last lines only: 



Pi = (pi/2,P3/2, ■ ■ -,^7+1/2, 2l_;^ ), Pn = { 2i_;^ ,Pn-J-1/2, 

j+i times j+i times 

Ui = {uq, Ml, ■ ■ ■ , uj, 0,--,0 ), Un-1 = ( p,-y)0 , Un-J, Un-J+I, 
j+i times j+i times 

Vectors 6aP and 6a^ are extracted from matrices 6D^p\6D^'^^: 

6aP = {6aQ, 6ai, 6a2, . . . , 6a^j, 6a^j, 6a^j_i, . . . , (5q;q)* 
6a^ = {6a^, fc", 5^2, • • • , Saj, 6a'}_^, . . . , ^Oq)* 



,PN -1/2) 

,un) (15) 



(16) 



It has to be noted, that operators P and U act from the space of the control variable a 
to the space of the model's solution u or p. Their matrices, consequently, are rectangular. 
Their dimensions are x 2( J + 1) and (A^ ~ 1) x 2( J + 1) respectively. 

So far, both 6a and {6u, 6p) are supposed to be small, we neglect their products in 
ffn]) and get 



d6u 

~dt 
d6p 

"dt 



(17) 
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with the same boundary conditions ([2]) for {5u, 5p). At initial time both 5u and 5p are 
taken to be zero because our study is confined at evolution of a pure perturbation due to 
boundary scheme. 

The same time stepping as in (Q is applied to (fT7|) : 



6u 



n+l 



6u 



n-1 



2t 



2t 



The first step of the tangent linear model ( fTTl) is written according to the scheme ( fTOl) . 
Taking into account the zero initial condition 5u{x, 0) = 0, 0) = we write 



2 ' ^ 2 

Equation (I18p can be rewritten in a matricial form: 



r , 



(19) 
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(20) 



with the first step (fT9|) 
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(21) 



To obtain the adjoint model for euclidean scalar product, we introduce adjoint variables 



and write backward evolution with transpose matrices (!20!) 
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The last step of the adjoint model is the adjoint of the first step of the tangent model: 




^(f/°)*(D(p))* 
r(pi/2)* 



/ 





l/2\* 




\ il } 



(24) 



where operators ([/")*, (D^^))*, (P")*, are adjoints to ^ and (HSl). 

We can see that the right hand side of the tangent linear model (|T71) is composed by 
two terms: or D^'^^Su and PSa^, or USa^. The first one, (fT2|) . is responsible for 

the evolution of a small perturbation by the model's dynamics, while the second one, 
(fTSll . determines the way how the uncertainty is introduced into the model. The first 
term is similar for any data assimilation, while the second one is specific to the particular 
variable under identification. This term is absent when the goal is to identify the initial 
point because the uncertainty is introduced only once, at the beginning of the model 
integration. But, when the uncertainty is presented in the approximation of derivatives 
near the boundary, or some other internal parameter of the model or of its numerical 
scheme, the perturbation is introduced at each time step. 

2.2. Cost function 

To perform variational data assimilation we introduce the following cost function: 

T 1 

u{a, X, t) - u°^%x, t)f + x, t) - p"^'{x, t)fdxdt = 



J(a) 





T 



u{a, X, t) — u°'"^{x, t) 
p{a, X, t) — p"^^{x, t) 



Mt 



(25) 



where the norm corresponds to Euchdean scalar product 



u{x, t) 
p{x,t) 



u{x, t) 
p{x,t) 



u{x, t) 
p{x,t) 



u'^{x,t) +p'^{x,t)dx (26) 



We suppose we have observations for all variables at any time. For numerical experiments 
in this paper we shall use the exact solution of the equation ([1]) as observations. This will 
help us to see the assimilation procedure and its results in the simplest and clear form. 
When this technique is applied to more complex model for which the exact solution is not 
available, we can use either real observations or the model's solution on a finer grid. 
To calculate the gradient of the cost function, we calculate first its variation 



61 = I{a + 6a) - I{a) = 

u{a, X, t) — u°^^{x, t) 
p(a, X, t) — p°''*(x, t) 







6u{x, t) 
6p{x, t) 



T 

2/«| 





M(a,x,t) - u°^'{x,t) 
p{a, X, t) — p°*'*(x, t) 




|» dt 



» dt 



8 



u{a, X, t) — u°^^{x, t) 
X, t) — p°^^{x, t) 



6aP 



|» dt 



(27) 



where T(t) ( 1 is the tangent model (j2Ti) . (j20l) integrated from t = to t and A{t) 



IS 



the adjoint model integrated from t to 0. 
Thus, the gradient of the cost function 



VJ 



A{t) 



u{a, X, t) — u"'^^{x, t) 
p{a, X, t) — p°^^{x, t) 



dt 



:28i 



is obtained as the sum of the adjoint model integrations. Each integration of the adjoint 
model starts from multiplication of the matrix ( l23i) by the vector 

/ u{a,x,i) -u°^'{x,i) \ 



p{a, X, t) — p°^'^{x, t) 





and followed by subsequent multiplications by matrices (123|) taken at corresponding time. 



This product is finally multiplied by the matrix 



to get the vector 



which 



represents the gradient of the cost function. 

This gradient is used in the minimization procedure that is implemented in order to 
find the minimum 

I{a) = minX(a) (29) 

Coefficients a are considered as coefficients realizing optimal discretization of the model's 
operators in the boundary regions. 

The minimization procedure used here was developed by Jean Charles Gilbert and 
Claude Lemarechal, INRIA 22|. The procedure uses the limited memory quasi-Newton 
method. 



3. Results of assimilation 

Exact solution of the equation ([1]) can easily be found by the method of variables 
separation. We look for solutions in a special form u{t,x) = a{t)b{x). A consequence is 
that 

— - — - -A 
a b 

The value of A is determined so that there exists a non-trivial solution of the boundary- 
value problem 

6"+A6 = 0, 6(0) = 6(1) =0 

Values of A are all positive, and the solutions are trigonometric functions. A solution 
that satisfies square-integrable initial conditions ([3]) for u and p can be obtained from 
expansion of these functions in the appropriate trigonometric series. 
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3.1. One trigonometric mode 

We shall analyze first the behavior of one trigonometric mode of the solution and 
further proceed with the analysis of more complex functions. 
Let us define the initial point for u and p in ([T]) as 

■u(a;, 0) = sin{knx) p{x,0) = cosiknx) (30) 

The solution of ([T]) determined by A = /c^vr^ is 

Uexact{x,t) = —\/2sm{k7it — tt/4) sin(A;7rx), Pexact{x,t) = -\/2cos(fc7rt — vr/4) cos(A;7ra;) 

(31) 

The solution fl3T|) will be used as artificial "observations" to be assimilated into the dis- 
cretized wave equation. The use of these data allows us to work with errors of numerical 
schemes only, avoiding all additional errors that may be present due to inexact data. 

Two numerical approximations are used for discretization of spatial derivatives in all 
internal points of the interval. Both discretizations are performed by formula ([5]), but one 
of them is of second order of accuracy with coefficients aj = (0,-1,1, 0) for j = (—1,0,1, 2) 
and the other one is of fourth order with Uj = ^(1, —27,27, —1). The simplest second 
order scheme on the boundary was used in both cases. That means both a" and in 
(j7]) were chosen to provide classical approximation of derivatives in points adjacent to 
boundary: 

'dp\ ^ P3/2 - Pl/2 ^^U\ ^ Ui - Up 

.dxJi h ' \dxJi/2 h 

In order to see precisions of these schemes we calculate the difference between the 
numerical solution u{x,t),p{x,t) and the exact one Uexactix,t),Pexact{x,t) and plot its 
norm 

1 

^{t) = J ({u{x,t) - Uexact{x,t)f + {p{x,t) - Pexact{x,t) f^dx. (33) 


Numerical solutions are obtained with k = 3, h = and r = j^o- 

It is well known that the principal error of classical (with approximations of derivatives 
near the boundary realized by fl32p ) solutions for both second and fourth order schemes 
consists in the wrong wave speed. Numerical solution of ([T]) is also composed of trigono- 
metric functions of the same amplitude but they oscillate with wrong frequency. The 
second order solution oscillates a little slower than the exact one, and the fourth order 
oscillates a little faster. 

In figiHA. and figUlB we see that the difference between exact and numerical solutions 
oscillates with the frequency Sir but have a growing amplitude. The velocity error is lower 
when the fourth order approximation is used, that's why the amplitude of the difference 
in figilB is lower than in fig{TJA.. 



10 



Figure lA. x — t diagram of the error of the Figure IB. x — t diagram of the error of the 

classical second order scheme. Contours from - classical fourth order scheme. Contours from - 
0.2 to 0.2 with interval 0.05. 0.04 to 0.04 with interval 0.01. 

If we look at figures figJ2l^ and 2B, we see the same phenomenon. The sohd hne 
in fig|2]A., that represents the norm of the difference between the exact solution and its 
second order numerical approximation, grows first up to value of 120 at time t = 108.3 
time units. After that, the norm decreases to at time t = 215.9 and restarts to grow. 
The fourth order approximation exhibits a similar behavior, the norm also grows up to 
value 120, but it reaches its maximum and the following zero at t = 491.1 and t = 982.2 
time units respectively. These moments of time, being beyond the picture window, are not 
shown. The speed error in the second order approximation results that at time t = 215.9 
numerical solution is exactly one wave period later than the exact one, and the difference 
between them vanishes. So far, the speed error is lower for the fourth order approximation, 
moments of the maximal and vanishing norm in figj2j3 are reached later. 




Figure 2A. Error £^{t) of the second order Figure 2B. Error ^{t) of the fourth order 

scheme: Classical - solid line, with assimilated scheme: Classical - solid line, with assimilated 
boundary - dashed line. boundary - dashed line. 

Thus, it was illustrated that the principal error of numerical approximation consists 
in the wrong wave speed. Indeed, if we apply numerical approximations to trigonometric 
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functions, we can calculate the error in the wave velocities. We substitute trigonometric 
solutions for u and p in the second order scheme, 

u{x,t) = sm{kx) sin{kt) = sin{ikh) sin{nkT) 
p{x,t) = — cos{kx) cos{kt) = — cos{ikh) cos^nkr) 



we get 



dpY _ P7+i/2-Pti/2 _ cos{{i + l/2)kh)-cos{{i-l/2)kh) 

COS( TZ/CT j " 



dxJi h h 

2 sin(^) sin{ikh) cos(zA;r) 
h 

/duy _ <+^-<-^ _^.^^.,^^^sin((n + l)A;r)-sin((n-l)A;r) 



sm(ikh)- 



\dtJi 2t ' ' 2r 

sin(/cr) sin{ikh) cos{ikT) 

T 

Thus, the first equation in ([1]) is approximated by 

duY fdpY /sin(A;r) 2sm{kh/2)\ 



sm{ikh) cos{ikT) 



hsm{kT) -2Tsm{kh/2) /dpY 

2rsin(A:V2) \d^)i ' 

Similar substitutions for u and p in the second equation give us the approximation of the 

system 



with 



<9p\" ^ /du\" , , 

-H^] =0 (35) 



dt ) i \dx 
dt J i \dx 

_ hsmjkT) 
- 2rsin(fcV2) ^ ^ 

Thus, we see that numerical wave velocity is equal to j32 rather than to one. 

If we perform similar manipulations with the fourth order spatial discretization, i.e. 
approximation of all spatial derivatives by (j5]) with stencil = ^(1,-27,27,-1), we 
get the velocity error 

12/isin(A;r) 

" 27rsin(A;/i/2) -rsin(3A;V2) ^ ' 

In figj3] we can see the form of speed errors — 1 and (34^ — 1 for three values of k. 
Horizontal axis is marked in values of ^. 
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Figure 3. Wave speed error (/? — 1) for second and fourth order 
schemes as functfons of r/h. 

We see that using second order scheme, we can simulate the exact solution. Indeed, 
when h = 2t the velocity of numerically approximated wave is exactly equal to the velocity 
of the theoretical solution for any wave-number k. Using any lower r we must assume the 
error in the waves velocity. 

On the other hand, it is impossible to calculate an exact solution with a fourth order 
scheme. The value of /34 — 1 vanishes in different points r/h for different k. The only 
conclusion we can make is the ratio r/h must either be small for this scheme, or some 
higher order time stepping should be used. 

For the given parameters {k = 3, h = and r = -^jj) errors in the wave velocity 
can be calculated by fl36|) and fl37|l : 



p2 = 3.09 X 10"^ Pi = -9.82 X 10^ 



These velocity errors determine the time when the numerical wave will be one period 
shifted with respect to the exact wave: T = ^^^^ Period _ ^j^^ second order 

scheme with k = 3 this time T is equal to 215.6 time units that corresponds well to 
numerically obtained 215.9. 

So, knowing errors produced by numerical schemes with chosen parameters, we shall 
perform the assimilation of the exact solution in order to see how these errors can be 
corrected by the optimal boundary discretization. 

We perform the data assimilation minimizing the cost function fl25l) assuming that the 
approximation of boundary derivatives is composed by two terms only (J in ([7]) is equal 
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to 1) and we get a numerical solution with no error in wave velocity. The norm fl33p of 
the difference between the exact solution and its optimal numerical approximation (lower 
lines in figj2j\ and B) oscillate around 3 x 10"'^ and 3 x 10~^ respectively. X — t plots of the 
difference u{x,t) — Uexact{x,t) presented in figJU^. and figlUB show very similar behavior 
of the error. The difference is composed of small moving waves that propagate back and 
forth between the boundary and the middle of the interval for both the second and the 
fourth order schemes. The amplitude of these waves is small comparing to errors of the 
classical scheme and, that is more important, remain small during any integration time. 
This fact can be seen in fig|2l Despite the data were assimilated during 6 time units only 
(T in ( l25i) is equal to 6), boundary approximation of derivatives has been sufficiently well 
identified to satisfy the model during any long integration, 300 time units and more. 




Figure 4A. x — t diagram of the error of the Figure 4B. x — t diagram of the error of the 

modified second order scheme. Contours from modified fourth order scheme. Contours from 
-0.03 to 0.03 with interval 0.01. -0.01 to 0.01 with interval 0.003. 

The choice of an optimal assimilation window (the time interval T during which the 
assimilation is performed) is obvious for this simple problem. Of course, T must not 
be too small. It must cover at least several wave periods in order to provide necessary 
information about errors in wave velocities. On the other hand, too long T is not optimal, 
because the assimilation over a longer interval is less efficient. First, we do not need too 
much data to assimilate because of the simplicity of the model. And second, too long T 
reduces computational efficiency of the method because of the necessity to run the model 
for a longer time in each iteration. 

Thus, assimilating the exact solution of the equation, we can construct an optimal 
approximation of boundary derivatives and obtain a rather accurate model which error is 
sufficiently small. However, boundary derivatives obtained in this procedure are strange 
from the point of view of approximation. 

When the second order approximation is used for derivatives in all internal points of 
the interval, the optimal discretization near the boundary obtainded by data assimilation 
has a form 



ox/ 1/2 



h ' 



dp 

dx 



3.014p3/2 - 2.828pi/2 
h 



{3t 



First of all, these formulas do not approximate a derivative. The first one approximates 
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the derivative multiplied by 1.048, the Taylor expansion of the second one has a form 



0.18^ + 2.92 

h 



(|)_,0.0m(g)_.0.I..^(g)^.O(.3) 



Neither expression for ^, nor for has arry reasonable order of approximation. 

^ OX OX 

The first one is of order, the second is of -1 order. Moreover, while we get always the 
same formula for approximation of the derivative of p varies in different assimilation 
experiments. Assimilations performed with different assimilation windows, for example, 

result in different coefficients for In fact, any combination , in (j7]) may be found 
as the result of assimilation under condition 

a? = -1.104ag- 0.107. (39) 

This linear relationship has been obtained experimentally performing assimilations 
with all assimilation windows in range from 600 to 2400 time steps (with the time step 
equal to 1/120 of the time unit). Resulting couples , presented in figj5]are positioned 
on a straight line with values varying from -1.5 to -5. 




Figure 5. Scatter diagram of ckq , obtained with different assimi- 
lation windows T in range from 5 to 20 time units. 

To explain these unusual approximations of the derivatives, we address first the u 
derivative, that is always approximated by 1.048 ^-^ . We know, the principal error of 
the classical scheme consists in wrong wave velocity. The data assimilation and control 
of the boundary derivatives can not modify numerical wave velocity. The only way for 
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this control to get a better solution consists in modifying the length of the interval. A 
numerical wave with wrong velocity will propagate on the interval with wrong length. 
But the length of the interval is adapted by data assimilation in order to ensure the wave 
with numerical velocity propagates the modified interval in the same time that the exact 
wave propagates the exact interval. So far, the control can not correct the error in the 
wave velocity, it commits another error in length in order to compensate the first one. 

As we have seen, the coefficient 1 in front of spatial derivatives in ([1]) has been replaced 
by (3 in fl35|) . Theoretical wave speed Cexact = 1 has, consequently, been replaced by 
numerical speed Cnum = P- The length of the interval Lexact = 1 should also be modified 
to satisfy 

Lexact LiYnodified j j ^exact Lexact 



modified J-'exact ~^ 
Cnum Cexact Cnum P 

However, the control can not modify all grid cells of the interval uniformly. It can act 
near boundaries only and can modify the length of cells just adjacent to boundary points. 
Hence, only two grid cells, one on the left and one on the right of the interval, can be 
modified. The modified interval, hence, becomes composed by A^ — 2 cells of length h = 
and two boundary cells of length 

2/.„o../.e. + (AT - 2)/. = ^ ^ = 1 - 1^ (41) 

For given parameters (N=30) = 1 + 3.09 x 10^'^ the boundary cells must be reduced to 
= (1 - 0.046)A. Consequently, the derivative m must be calculated over 
modified cell 



du\ ui-uq 1 ui-uq 1 048^^""° 



.dxJi/2 hmodified 1 - 0.046 h h 
This is exactly the coefficient obtained in the data assimilation for the derivative of u 



So, we can state that it is reasonable to obtain wrong approximation of derivatives 
near boundaries as a result of data assimilation. This error compensates the error of the 
wave speed. 

As for derivatives of p, they must also be modified. The only difference with u consists 

in fact that is calculated over two half of cells: one half of the first cell (adjacent 

to boundary point), and one half of the second one, next to the first. Hence, only one 
half of the modified cell participates in the derivative of p and its modification is 

dp\ _ P3/2-P1/2 _P3/2-pi/2 2/i ^^2) 



dxJl hmodifiedf^ + h/2 h hmodified + h 

In this experiment we should have obtained = 1.023 ^'^^^ ^ ^^^^ . 

And indeed, the couple = —1.023 , = 1.023 belongs to the set (l39l) . This is the 
only point on this line where + = and the derivative is approximated with zero 
order rather than minus first order. 

Non uniqueness of optimal and ag can be explained if we take into account that p 
has also a form of cosine of Snx. Hence, at any time pi/2 = A{t) cos(37r/;,/2) and P3/2 = 
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A{t) cos(97r/i/2) with some A depending on time. Their hnear combination aiPi/2 + (yoP3/2 
can vanish if 

p 

" ~4cos2(A;7rV2)-3" ^^^^ 

Consequently, all couples , belonging to the line that passes by the point = 

— 1.023 , a? = 1.023 with tangent —-. ^, = —1.108 produce the same deriva- 

' ^ ^ 4cos^(37r/i/2) - 3 ^ 

tive. This line coincides withing accuracy of computation with the set (139!) obtained 
numerically. Any point on this line gives coefficients ap that theoretically provide the 
same value of the derivative and the same value of the cost function. This line forms the 
kernel of the Hessian of the cost function. 

Numerical approximation of the solution is slightly different from cosine and numerical 
approximations of the derivative obtained with different coefficients from the kernel 
are not exactly the same. The assimilation chooses the best fitting point in the kernel for 
particular experiment that provides slightly lower value of the cost function. The choice 
of this point depends on particular parameters of the experiment such as assimilation 
window. That's why we get different pairs ag, in different experiments. All these 
pairs are in the kernel of the Hessian, they provide almost the same cost function values, 
but each of them corresponds better to one particular window. If we are interested in 
optimal boundary scheme for the whole model rather than in the best fitting point for 
a given assimilation window, we may define another criterium of choice and impose this 
criterium in the cost function. One choice, usually assumed in data assimilation, requires 
that optimal point must be situated not far from the initial guess. However, adding this 
requirement would not allow us to choose one point in the kernel. The requirement of 
low distance from the start would draw the optimal point out of the kernel because, as 
we have seen above, the initial guess point is not situated in the kernel. 

Instead of imposing low distance from the starting point of minimization, we prefer to 

require the term in the Taylor expansion with the order minus one to be equal to zero. 

J 

This implies the sum ctj = must vanish. For this purpose we add the term 

i=o 

R = v{j:a,y (44) 

j=0 

to the cost function fl2Sl) and appropriately modify its gradient 0281] adding the term 

J 

VR = 2r]Y,aj. (45) 

j=0 

Imposing sufficiently large weight 1] we get the only approximation of p derivative for any 
assimilation window. The derivative is approximated by ag = —1.023 ,ai = 1.023 that 
ensures vanishing first term in the Taylor development. 

Modification of the cost function by has a very small infiuence on the final value 
of the cost function because this modification determines the choice of the particular point 
in the kernel of the Hessian. 

Finally, we note that there is no significant difference in the final value of the cost 
function in experiments with different J in ([7]). Several experiments have been carried out 
with 2, 3 and 5 controlled coefficients a, but the minimization procedure has converged 
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always to the same value. Obviously, two control coefficients are already sufficient in 
this simplest case. Adding supplementary a just increases the kernel dimension with no 
influence on the cost function. 

3.2. Two trigonometric modes 

When initial conditions of the model ([T]) are more complex than one trigonometric 
mode, the exact solution of the wave equation is a linear superposition of exact solutions 
corresponding to each trigonometric mode of the Fourier development of initial conditions. 
Each mode has it's own frequency and propagates with it's own velocity. 

Numerical solution for each Fourier mode commits an error in the wave velocity. But, 
as it has been discussed above, this error is different for different modes because it depends 
explicitly on the wavenumber k (l36l) . (!37|) . Consequently, in presence of multiple Fourier 
modes, the interval length must be modifled in order to correct different errors in wave's 
velocities simultaneously. 

We consider flrst a superposition of just two waves with k = 2 and k = 5. We see from 
the equation fl4ip that to compensate the error in the wave velocity for the wave with 

k = 27r, the control must modify the length of the boundary cell by '^"'j^f'""^ = i _ 0.020 
and the coefficient in front of the approximation of the derivative of u at point 1/2 must 
be 1.021. In the same time, the velocity error for the wave with = Svr is compensated 
^j^gn hmodified ^ ^ _ g^^28 and the coefficient in front of the derivative 1.142. 

Performing experiments with both wavenumbers k = 2 and k = 5 separately and 
with their superposition, we see in flgEJ\ that the data assimilation procedure is able 
to compensate the error in wave velocity in all three cases. The cost function of the 
model with original coefficients shows wrong velocities of numerical waves in all three 
experiments, but the model's solution with optimal coefficients is much closer to the 
exact one. We see the cost function values as low as 3 x 10~^ for the wave with k = 2 
and 10~^ for the wave with k = 5. The line that corresponds to the cost function in the 
experiment with two waves superposed is indistinguishable from the line corresponding 
to the experiment with k = 5. They oscillate both around X = 10^^. That means the 
residual error of assimilation of the superposition of two waves is close to the biggest error 
of assimilation of each particular wave. 

In order to analyze the expression that is used to calculate the derivative of u near 
boundaries in flgl6]B, we perform a set of assimilations with all assimilation windows in 
range from 600 to 2400 time steps (with the time step equal to 1/120 of the time unit) 
for all three types of initial conditions of the model, i.e. one wave with either k = 2tt 
or k = and both of them. When /c = 27r we get always the same resulting couples 

= —1.021 , a" = 1-021 as expected. Coefficients , in the experiment with k = 5n 
are also all positioned near the theoretical value ±1.142, but not as concentrated as in the 
experiment with k = 2tt. Values in this experiment are distributed in the interval from 
1.138 to 1.144. Obviously, the wave length of the wave with A; = Svr is too short to be well 
reproduced by a 30 points resolution grid. This coarse resolution adds numerical noise in 
the solution and leads to the dependence of the assimilation result on the window. 

Optimal coefficients , a" in the experiment with two waves are situated in the 
middle of the flgure flgl6]B. We can note two particularities. First, their distribution is 
even more dispersive than with k = Svr: they occupy the interval from 1.07 to 1.09. And 
second, expressions for u derivative near the left and near the right boundary are no 
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longer the same. One can see in figEP the set of coefficients , a" in this experiment is 
sphtted into two subsets with a gap between them. 



DiH'erence with the exact solution 
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Figure 6 A. Cost function of numerical so- 
lutions for modes with k = 27r, k = 5tt and 
their superposition. Solutions with the original 
boundary scheme are plotted with solid lines, 
with identified schemes - by dashed lines. 



du 



Figure 6B. Optimal coefficients for ^ ap- 
proximation for modes with k = 2-k (lower right 
corner), fc = 57r (upper left corner) and their 
superposition (center). 

Coefficients of expressions for p derivative in the experiments with two waves (not 
shown) possess also a kernel that form the line situated between lines obtained in exper- 
iments with single waves. 

3.3. Other functions 

If we consider an arbitrary functions as initial conditions of the wave equation, we 
have all admissible Fourier modes in the solution. In order to see the action of the control 
in this situation we perform the data assimilation for the model with initial conditions 
prescribed as 

u{x, 0) = 20x^(1 - x)e-^^', 0) = (x - 0.5)6^"^ (46) 

Combining polynomials and exponents we ensure that different trigonometric modes are 
present in the spectrum of initial data that leads to a rich spectrum in time. 

First of all, the control of just two coefficients in expressions for derivatives is no 
longer able to ensure non growing cost function beyond the assimilation window. We see 
in fig|7|A. that the cost function of the model with optimal coefficients , ai grows after 
the assimilation end in the same way as the cost function of the original model. Solid and 
upper dashed (that corresponds to J = 1) lines are parallel to each other. In fact, the 
data assimilation reduces the model's error approximately 20 times, but the behavior of 
the error remains the same. Consequently, we can not state that the model's error with 
optimal boundary approximation will always be small. Increasing with time, the error 
will later reach the same values as the error of the original model. 
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Difference with the reference solution 
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Figure 7A. Cost function of numerical so- 
lutions with J = 1 and J = 4. Solutions with 
the original boundary scheme are plotted with 
solid lines, with optimal scheme - with dashed 
lines. 

This fact can be explained by the analysis of the expression PX 
scheme 



Figure 7B. First two optimal coefficients 
for ^ with J = 4 in experiments with different 
assimilation windows. 



for the second order 



''modified 



h- 



Nhp-1 
~2 ~ 



h- 



2/3 



h 



r sm{kh/2)\ 
h sin(A;r) / 



The coefficient in the expression for the derivative of u becomes 



h 



h? sin(A;r) 



hmodified {h - h/2) sin(A;r) + r sm{kh/2) 



(47) 



(4J 



For the given parameters (/i = 1/30, r = 1/120) we get c = 3^5 cos(fc/120) 14 

denominator of this expression vanishes and changes sign when k 



The 



14.0267r. Consequently, optimal expression for ^ at the first point for the wave with 
k = ISvr must have an opposite sign with respect to the classical approximation, namely: 



120arccos(||) 



du 
'dx 



1/2 



-7.05^^1^. 



The wave with k = ISvr is present in the spectrum of initial 

conditions (its wavelength is equal to ^ = 4/i) but corresponding optimal expression for 
the derivative can not be obtained in the assimilation procedure because the scheme is 
instable with negative c. Hence, the minimum is unreachable and we can not obtain the 
optimal approximations of derivatives near the boundary. Data assimilation allows us to 
compensate the error in wave velocities for first 14 trigonometric modes, but all other 
modes continue to propagate with wrong velocities. That's why the cost function in the 
experiment with assimilated data is smaller than the original cost function, but the long 
time behavior is similar in both experiments. 

In order to obtain the cost function that does not increase after the end of assimilation, 
we may try to control more coefficients a in ([71) in order to be able to identify optimal 
coefficients in the domain where the scheme is stable. Increasing the number of controlled 



20 



parameters, we increase the number of degrees of freedom and the dimension of the kernel 
of Hessian. The intersection of the kernel and the region where the scheme is stable may 
become non null and allow the assimilation to reach the minimum. 

Indeed, if we perform assimilation with J = 4, i.e. 5 coefficients a in ([7j), we get 
smaller non increasing cost function (lower dashed line in fig|7|A.). 

Coefficients a" in the experiment with J = 4 are distributed in a wide area, showing 
larger multidimensional kernel of the Hessian. An example of such a distribution is shown 
in figI7]B. To obtain this figure, we perform a set of experiments with different assimilation 
windows in range from 800 to 5000 time steps of the model. In each assimilation we get 
different sets of coefficients a but almost the same final cost function showing all obtained 
a are in the kernel of the Hessian. Only the first two coefficients are plotted in figI7]B- One 
can see, they occupy much wider area than in experiments with one or two trigonometric 
waves and J = 1 shown in fig|6j3. 

4. Conclusion 

The purpose of this paper is to study the variational data assimilation procedure ap- 
plied for identification of the optimal parametrization of the derivatives near the boundary 
on the example of a simple wave equation in view to use this kind of data assimilation in 
ocean models. Consequently, conclusions are formulated from this point of view. 

Comparing this procedure with now well developed data assimilation intended to iden- 
tify optimal initial data, we can say there are both common points and differences as well. 

Tangent ffT7|) and adjoint fl25]) models are composed by two terms, presented by f lT^ 
and (fTSll . The first one, governs the evolution of a small perturbation by the 

model's dynamics. This term is common for any data assimilation no matter what pa- 
rameter we want to identify. The second one, U or P, (ITSl) . determines the way how the 
uncertainty is introduced into the model. So, if we intend to identify an optimal bound- 
ary parametrization for a model with an existing adjoint developed for data assimilation 
and identification of initial point, we can use this adjoint as ( I12p part because this part is 
common for any data assimilation. However, the part decribed by f|T5l) must be developed 
from the beginning because it is specific to the particular control parameter. This develop- 
ment may be technically difficult for complex models, especially on grids with distributed 
variables like Arakawa's "C"-grid. Numerous interpolation and differentiation operators 
are frequently applied successively to a model's variable on these grids resulting in non- 
linear dependence of the model's state on control coefficients. Development of the adjoint 
model and, particularly, it's (fT5l) part, is complicated by working with nonlinearities of 
higher degree. 

Another difference consists in the number of control parameters and their dimensions. 
The dimension of initial point of the model is usually equal to the dimension of the 
model's state variable. Contrary to this, when we control boundary parametrization, the 
dimension of control variables is very different from the dimension of the model's variable. 
Moreover, the dimension of the control might be lower than the dimension of the model 
state because the dimension of the control is proportional to the length of the boundary 
of the domain, while the dimension of the model's state relates to the area of the domain. 
That means the quantity of controlled parameters and the dimension of the gradient 
of cost function may be much lower than the quantity of variables in the models state. 
Taking into account mentioned technical difficulties in development of the adjoint, it may 
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be reasonable to try to calculate the gradient by some other method beginning with the 
simplest finite difference method. Of course, this will be more expensive computationally, 
but the gain in the development procedure may compensate this excessive computational 
cost. 

Concerning the data assimilation results, we see the data assimilation can correct 
errors of numerical scheme by controlling approximations near boundaries. This fact 
may be very useful in applications of this method to the ocean models. In addition to 
natural corrections of the position of the rigid boundary and prescribed physical boundary 
conditions, we may hope to be also able to improve the quality of the scheme that is used 
in internal points. 

We can see in these assimilation experiments the presence of a kernel of the Hessian. 
Consequently, the choice of optimal boundary parametrization is not unique. However, 
all sets of coefficients a from the kernel are equivalent: they provide the same (or almost 
the same) cost function's value and almost the same evolution of the model's solution 
after the end of assimilation. In the same time, we can note that optimal parametrization 
of derivatives near the boundary may approximate nothing in classical sense, i.e. it may 

not be valid for an arbitrary function. We have seen here that obtained expression for ^ 
is valid for the cosine-type functions with appropriate wavelength only. Hence, we must 
take into account that coefficients found by data assimilation are valid for given model's 
parameters only. 

In the last experiment in this paper, with the wave composed by multiple trigonometric 
modes, we have encountered the necessity to increase the number of control parameters. 
In the case when the optimum is unreachable, increasing the kernel dimension allows to 
obtain better results. Combining the number of controlled coefficients (that increases the 
kernel dimension) and the possibility to dump the first term of the Taylor development 
of the resulting expression by (jHj) (that decreases the kernel dimension) may help us to 
get a reasonable result. 
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